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^ ' Abstract 

o. 

Qs , We derive an analytical trace formula for the level density of the two-dimensional elliptic billiard 

CTs ' using an improved stationary phase method. The result is a continuous function of the deformation 

parameter (eccentricity) through all bifurcation points of the short diameter orbit and its repetitions, 
and possesses the correct limit of the circular billiard at zero eccentricity. Away from the circular 
limit and the bifurcations, it reduces to the usual (extended) Gutzwiller trace formula which for 
^ ' the leading-order families of periodic orbits is identical to the result of Berry and Tabor. We show 

^ , that the circular disk limit of the diameter-orbit contribution is also reached through contributions 

G ' from closed (periodic and non-periodic) orbits of hyperbolic type with an even number of reflections 

from the boundary. We obtain the Maslov indices depending on deformation and energy in terms 
of the phases of the complex error and Airy functions. We find enhancement of the amplitudes 
p^ ' near the common bifurcation points of both short-diameter and hyperbolic orbits. The calculated 

J~j , semiclassical level densities and shell energies are in good agreement with the quantum mechanical 

ones. 
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1 Introduction 

The periodic orbit theory (POT), developed by Gutzwiller [1, 2] for chaotic systems, by Bahan and 
Bloch [3] for cavities, and by Berry and Tabor [4, 5] for integrable systems, has proved to be an important 
semiclassical tool not only for an approximate quantization, but also for the description of gross-shell 
effects in finite fermion systems [6, 7]. Gutzwiller 's approach has been extended to take into account 
continuous symmetries [6, 8, 9, 10, 11, 12] and is therefore applicable to systems with mixed classical 
dynamics, including the integrable and hard-chaos limits. 

An important role is played by the classical degeneracy of the periodic orbits in systems with continu- 
ous spatial or dynamical symmetries: the orbits are then not isolated in phase space (as it was assumed in 
Gutzwiller's original trace formula, and as is the case in chaotic systems), but occur in degenerate families 
with identical actions. The degree of degeneracy /C is defined as the number of independent parameters 
which are necessary to uniquely specify an orbit within each family. E.g., the orbit families with the 
highest degeneracy in spherical systems with spatial 5*0(3) symmetry have /C = 3, corresponding to the 
three Euler angles that specify the orientation of an orbit within the plane of motion and the orientation 
of the plane itself; the orbit families in two-dimensional systems with U(l) rotational symmetry have 
JC = 1: the isotropic harmonic oscillator in 2 dimensions has SU{2) symmetry and hence orbit families 
with JC ~ 2. Orbits with different degeneracies /C may also occur in one and the same system, such as 
the spherical cavity discussed by Balian and Bloch [3] where the diameter orbit has /C = 2 and all other 
orbits have /C = 3; the spheroidal cavity [13] where /C = 2, 1 or occurs (the latter corresponding to 
isolated orbits); or the elliptic billiard with /C = 1 or 0, as discussed in the present paper. 

However, problems arise for all these trace formulae in connection with the breaking of a continuous 
symmetry and with the bifurcation of stable periodic orbits when a continuous parameter (energy, defor- 
mation, external field) is varied. The reason is that at such critical points the standard stationary phase 
approximation, used for integrations in the derivation of the trace formula, breaks down and leads to 
divergences and/or discontinuities of the amplitudes in the trace formula. This happens most frequently 
in mixed systems, but it occurs also in integrable systems. Typical examples are the two-dimensional 
elliptic billiard and the three-dimensional spheroidal cavity. In the former, all repetitions of the short di- 
ameter orbits undergo bifurcations at specific deformations, whereby new families of hyperbolic orbits are 
created. Similarly in the latter system, the periodic orbits lying in the equatorial plane perpendicular to 
the symmetry axis bifurcate also at specific deformations, whereby new 3-dimensional orbits appear [13]. 
In both systems, all bifurcations and the limit to the spherical shape lead to divergent amplitudes in the 
trace formulae, see Refs. [6, 11, 14, 15, 16, 17, 18, 19, 20, 21]. Since for each family with a given value 
of /C, the extended Gutzwiller trace formula [6, 8, 9, 10] has an amplitude proportional to h^^^^ ''^' , it 
is evident that the breaking of a continuous symmetry must be accompanied by a discontinuous change 
of the amplitudes, which manifests itself in the form of a singularity when one attempts to reach the 
unbroken symmetry limit. (An exceptional situation occurs in anisotropic harmonic oscillators when 
changing from irrational to rational frequency ratios: here the divergences of the different periodic or- 
bit contributions have been shown [23] to cancel identically, such that the trace formulae — which are 
quantum-mechanically exact here — hold for arbitrary frequency ratios, although their analytical form 
changes in the different limits; see also Rcf. [7].) 

Since symmetry breaking and orbit bifurcations occur in almost all realistic physical systems, there is 
a definite need to overcome these singularities. The importance of bifurcation effects in connection with 
the emergence of the 'superdeformed' shell structure in atomic nuclei was emphasized in Refs. [6, 18, 20, 
21, 22]. In order to improve the POT in these critical situations, various methods have been proposed. 
Like in the treatment of continuous symmetries considered in Refs. [8, 9, 10, 11], they essentially consist in 
taking some integrals in the derivation of the trace formula more exactly than by the standard stationary 
phase method (SPM). 

Berry and Tabor suggested in Ref. [4] a quite general method to treat bifurcations in integrable 
systems. Starting from the trace integral for the level density in action-angle variables, they reduce it 
to the Poisson-sum trace formula and do all trace integrations except one by the SPM, extending the 
integration limits from — oo to +oo. At bifurcations, this leads to singularities in the amplitudes when the 
stationary points are close to the limits of the integration range. According to Ref. [4], in this case one 
has to take the integral within the exact finite range. The integration range need not necessarily include 
the stationary points (in the case of negative or complex stationary points), but the latter are assumed 
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to be close to the integration limits. For intcgrable systems, this idea was applied to the periodic-orbit 
families with the highest degeneracies, for which one can exactly carry out the integrals over the action 
angles, giving 27r for each degree of freedom [5]. This was the starting point of a uniform approximation 
that was further developed by various authors [24, 25, 26]. 

Another type of uniform approximation was initiated by Ozorio de Almeida and Hannay [27] (see 
also Ref. [28]) and further developed by Sieber and Schomerus [29, 30, 31] for various generic types of 
bifurcations. Writing the trace integral in a phase-space representation, they expand the action around 
the bifurcation points into so-called normal forms which usually can be integrated analytically with finite 
results. The correct asymptotic recovery of the Gutzwiller amplitudes far from the bifurcation points 
can be obtained by a suitable mapping transformation whereby the amplitude function, together with 
the Jacobian of the mapping transformation, is expanded up to an order consistent with that of the 
action in the exponent of the integrand. Near the bifurcation points, there is a common contribution of 
all participating (real or complex, so-called 'ghost') orbits to the trace formula. 

A similar technique, starting from the Berry- Tabor approach for intcgrable systems and using a 
'pendulum mapping', was used by Tomsovic, Grinberg and UUmo [32, 33] to derive a generic uniform 
approximation for the breaking of orbit families with a one-dimensional degeneracy, corresponding to 
U{1) symmetry, into pairs of stable and unstable isolated orbits. Finally, some analytical uniform trace 
formulae for the breaking of the higher-dimensional SU{2) and 50(3) symmetries in specific two- and 
three-dimensional systems have been derived very recently [34] . Hereby the trace integral was performed 
over the de Haar measure of the corresponding symmetry groups, as in the derivation of the unperturbed 
trace formulae for these continuous symmetries, [10] and the mapping was done onto the forms of the 
action integrals obtained in perturbation theory [35, 36]. 

It should be mentioned that all the uniform approximations mentioned above can be used only for one 
isolated critical point of symmetry breaking or orbit bifurcation: they fail, in particular, [29, 30, 31, 33, 34] 
when two critical points arc so close that the actions of the participating orbits at these points differ 
by less than ~ h. To our knowledge, no common uniform treatment of two nearby bifurcations (in the 
above sense), or of a bifurcation near a symmetry-breaking point, has been reported so far. 

In this paper, we propose an approach to simultaneously overcome the divergences due to symmetry 
breaking and any number of bifurcations in the two-dimensional elliptic billiard and the three-dimensional 
spheroidal cavity. Although our framework is quite general, we limit here its application to the elliptic 
billiard. The three-dimensional spheroidal cavity will be treated in a succeeding paper, [13] and the 
extension to non-integrable systems is planned for future research. We start from a phase-space trace 
formula, [11, 37] which after some transformations becomes identical to that obtained from the mixed 
phase-space representation of the Green function in Refs. [30, 38], as explained there and further below in 
this paper (see §4.3). Analogous versions of the phase-space trace formulae were suggested in Refs. [5, 10]. 

In contrast to previous investigations, [4, 5, 24, 25, 26] we calculate the integrals over angles, too, 
by the stationary phase method. Note that we also include orbits with lower degeneracies, such as 
the isolated diameters in the elliptic billiard and the equatorial orbits in the spheroidal cavity, hereby 
extending the method of Ref. [4] . Our main point is that the stationary-phase integrals over both action 
and angle variables are calculated with expansions of the phase and amplitudes like in the standard SPM, 
but within finite intervals in all cases where it would lead to divergences if one or both integration limits 
were taken to oo or — oo. We will also discuss the role of non-periodic closed orbits (see §5.4). For the 
Maslov indices, which for the bifurcating orbits depend on the deformation and near the critical points 
also on the energy, we follow the basic ideas of Maslov and Fedoryuk [39, 40, 41, 42]. We obtain separate 
contributions to the trace formula from the bifurcating periodic orbits, and we remove the singularity of 
the isolated long diameter (i.e., the separatrix) near the circular shape of the elliptic billiard in a simpler 
way than in Ref. [26] . 

In this way we obtain an analytical trace formula for the elliptic billiard which gives finite and con- 
tinuous contributions at all deformations, including the circular disk limit and all bifurcation points of 
the short diameter orbit. Although its derivation and its explicit form are quite different, our final trace 
formula is similar to the uniform approximations mentioned above in the sense that it connects smoothly 
to the standard (extended) Gutzwiller trace formulae for the different orbit types for deformations suf- 
ficiently far away from all critical points. 
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2 Phase-Space Trace Formula in the Closed Orbit Theory 

2.1 Semiclassical trace formula 

The level density g(e) is obtained from the Green function G{r' ^r" ]e) by taking the imaginary part of 
its trace: 



g{e) ^--Im f dr" f dr' Gir' , r"; e)6{r" - r') 
= --Im f dr" f dr' f dpGir' ,r";e) exp 



{p-{r"-r') 



(1) 



Within the semiclassical Gutzwillcr theory, [1, 2] the Green function G(r' ,r";e) can be represented in 
terms of the sum over all classical trajectories a connecting two spatial points r' and r" at fixed energy 
e. Inserting it into (1), we obtain the semiclassical level density 

g^ci(g) - (2,;,)(L+i)/2 lmY^Jdr"JdpJdr'\J{p\t^y\er' 

xcxp{^[^„(r',r",e)-p.(r"-r')]-yA^o}. (2) 

Here Sa{r' ,r" ,e) ~ J , dr ■ p is the action along the trajectory a, n is spatial dimension, and /Iq, is 
related to the number of conjugate points (i.e., turning and caustics points along the trajectory) [42]. 
J^aip' ,ta',r" ,s) is the Jacobian for the transformation from initial momentum p' (at the point r') and 
time interval ta (for the classical motion along the trajectory from initial to final point) to final coordinate 
r" and energy s. 

2.2 Phase space variables 

Integrating over r' in Eq. (2) along the direction transverse to the trajectory a by the stationary phase 
method (SPM), we are left with the integral over the component of dr' parallel to the trajectory, which 
gives just an energy conserving delta function d{e — H{r',p')). We hence arrive at the phase-space trace 
formula [37] 

X exp{^ [So.{p',p",t^) + [p" -p') ■ r"] - iv^y (3) 

Here i/(p" ,p^) is the Jacobian for the transformation from initial to final momentum components p'^ 
and p" , respectively, perpendicular to the trajectory a. This Jacobian is equal to one of the elements of 
the stability matrix (see, e.g., Ref. [7]). Sa{p' ^p" ,ta) is the action in the momentum representation 

5a(p',p",ia) = - r dp -rip), (4) 

Jp' 

which is related to the usual action in coordinate space 

5„(r',r",£)= / dr-p{r) (5) 

J r' 

by the Legendre transformation 

5„(r',r",e) -p' • (r" - r') = S^{p',p",t^) + [p" - p') ■ r" . (6) 

Note that the integrand in the phase-space trace formula (3) (except for the exponent related to the 
phase part proportional to r") is the semiclassical Green function in the mixed representation which 
contains explicitly an energy-conserving (5- function in our case, unlike the form discussed in Ref. [10]. 
(Consequently, the momentum components are not independent, which is important for the following 
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application of the stationary phase method; see more details in the next subsection and in §4.) Due to 
energy conservation, i.e., H{r' ,p') = H{r" ,p"), the trace formula (3) can be rewritten in an alternative 
form where the integration variables are changed from {r" ,p') to [r' ,p"). The sum in (3) runs over all 
isolated classical trajectories a with starting momentum p' and final point r" (or with starting point 
r' and final momentum p" in the alternative form) , for a fixed time interval t^ of the classical motion 
along a. 

2.3 Periodic orbit theory 

The trajectories a in the phase space trace formula (3) arc not necessarily closed orbits in the usual 
coordinate space. But after separation of the extended Thomas-Fcrmi part (corresponding to the 'zero 
length orbits') and integration over one of the momentum components exploiting the (5-function, we 
shall use further scmiclassical approximations. We first write the stationary-phase conditions for the 
integration variables in (3). The stationary conditions for the momentum variable p' are the closing 
condition for the trajectories a in the usual coordinate space, r' — r", and the Jacobian in Eq. (3) is 
unity due to the Liouville theorem of the phase-space volume conservation, see Rcf. [7] . The additional 
stationary-phase conditions for the integration over spatial variables r" selects the periodic orbits, p' = 
p" , and we obtain the POT and all known trace formulas including the Poisson-sum trace formula [37]. 
We shall then integrate over components of the phase-space variables exactly if we have identities for 
them. Other integrations will be done by an improved stationary phase method (ISPM). Tmprovcd' here 
means that we carry out the integrations in finite ranges, after expanding the exponent of the integrand 
around the stationary point up to second order terms, and taking the amplitude at the stationary point 
(or use a higher-order expansion of amplitude and phase, if necessary). All stationary points which appear 
outside the physical region of the integration over the phase-space variables are also taken into account, 
even if they are complex. In this way we get simple and continuous analytical solutions that stay finite 
at all critical (bifurcation and symmetry-breaking) points. Different from other uniform approximations 
mentioned in the introduction, our results appear as explicit sums over separate contributions that 
correspond to the periodic orbits in the asymptotic regions away from the critical points. 

3 Classical Mechanics 

3.1 Elliptic billiard as integrable system 

We consider an elliptic billiard with axes a and b (with a < b) along the y and x coordinate axes, 
respectively, and ideally reflecting walls. This is an integrable system which can be separated in the 
elliptic coordinates {u,v) defined in terms of the cartesian coordinates {x,y) by 



a; = Ccosusinhu, y — ^sinucoshw, C~\/b'^ — a'^, (7) 

with 

-| < « < |, 0<v<Vb. (8) 

Hereby {x,y) — (±C, 0) are the foci of ellipses given by w = const., and v — Vb is the elliptic boundary. 
It is convenient to introduce the deformation parameter 77 — b/a > 1 and to keep the area of the ellipse 
constant by setting ab — R^, so that one gets b = R^ and a ~ R/\/V- The second constant of the 
motion, besides the energy e, is the product of the angular momenta I- and 1+ with respect to the two 
foci. For the following, it is advantageous to use the single- valued quantity a defined by 

There are two types of orbits, depending on the relative sign of Z_ and ly. elliptic orbits circulating 
around both foci for l_l^ > or cr > 1, and librating hyperbolic orbits for l_l^ < or cr < 1. Their 
names used here indicate that the former are limited to the area between the elliptic boundary given by 
V = Vh and a confocal elliptic caustic given hy v ~ Vc, whereas the latter are confined to the area between 
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Figure 1: Some classical periodic orbits in the elliptic billiard are shown by thin solid lines. Left-hand 
side; elliptic triangular (1,3) and rhomboidal (1,4) orbits: Right-hand side; hyperbolic butterfly orbit 
(1,4), fromRef. [11]. 



the two branches of a hyperbolic caustic given by u = ±Wc and the elliptic boundary. The critical values 
for the boundary and the caustics are given by 



Vb = v/vV^ ^ Ij Vc = arccosh(l/\/a'), Uc — arcsin(\/CT). 
In terms of the above quantities, the single- valued action integrals /„ and ly become 

pC 



(10) 



Pudu - 



sin u, 



du y a 
ly = (pPvdv = — / dv Vcosh V ~ a, 

J 71" Jv, 



— Uc 



(11) 



where p — \/2me — hk is the constant classical momentum of the particle. Since the system is integrable, 
its Hamiltonian depends only on the actions and not on the variables u, v, i.e., H{Iu, ly, u, v) = H{Iu, ly)- 



3.2 Periodic orbits 

As shown by Berry and Tabor, [4] the periodic orbits of an integrable system are found by the condition 
that the angular frequencies (for angle variables conjugate to the actions) have rational ratios. In the 
present case, these frequencies are given by tOu = dH/dlu, i^v — dH/dIv, so that the periodic orbits are 
characterized by pairs of positive integers (M„, My) 



UJy 

UJv 



1- 



F(g,^) 
F(f,«) 



where 



K == sin Uc/ cosh Vc 



^, (M„ > 1, M, > 2M„), 

IVlv 



arcsin(cosh Vc/ cosh Vb ) , 



(12) 



(13) 



and F{d,x) is the elliptic integral of the first kind [47]. The greatest common divisor of M„ and My 
corresponds to the repetition number M = 1, 2, 3, . . . of a primitive periodic orbit (n„, riy): 



{Mu,My) = {Mnu.Mny) = M{nu,ny). 



(14) 



The solutions of Eq. (12) for k and 9 which correspond to families of degenerate periodic orbits with 
/C = 1 are, labeled accordingly for elliptic and hyperbolic orbits. 



i^h 



1 



arcsm 



3in (Va(l - 1^2)) 



(15) 



Figure 1 shows the shortest periodic orbits of each kind. The degeneracy parameter K, was defined as 
the number of parameters that specify the orbits within a family with a common action. Due to the 
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separation of variables in elliptic coordinates (7) we have two single-valued action integrals /„ and ly 
(11). They are related through the energy conserving equation e = H{Iu,Iv), and can be written in 
terms of one parameter of the family a (or 1-1+), i.e., we have /C = 1 (see Refs. [6, 8, 9, 11, 50] for more 

details.) 

3.3 Energy surface 

For the energy surface e — H{Iu,Iy) one can get from Eqs. (11) the parametric equations (83) for the 
elliptic orbits and (84) for the hyperbolic orbits [19]. The energy curve (83) or (84) can also be considered 
through the single-valued parameter a or double-valued k defined within the same range < k < 1 for 
both kinds of orbits. The solutions a found from the periodic orbit equations (12) for elliptic orbits 
satisfy the inequality cr > 1 in the elliptic part (83) of the energy curve. On the other hand, cr < 1 
for the hyperbolic part (see Fig. 2a). The two regions are separated by the separatrix point as — 1, 

(s) 

corresponding to the long diameter orbit, where the value of the action I^ — lu is given by 

/(^) = 2K/7r. (a, === 1) (16) 

Thus, each phase space torus is split into two regions by the separatrix: a hyperbolic and an elliptic 
region. In the hyperbolic part (0 < cr < 1), the action variable /„ changes from to the separatrix value 
/„ . In the elliptic part (1 < cr < cr^r), lu changes from the separatrix value to the maximum value lu 
that corresponds to a 'creeping' (or 'whispering gallery') orbit and is given by 



j(cr)_ 2pi^V^ ^/^^ 1 \_ 2pR^ I^ ,/^f^l 



TT \ 2 ^CTcr / TT \ 2, V I 

(j„ = cosh^ Vb = ?7^/(?7^ - 1). (17) 

The short diameter (1,2) and its repetitions M(l,2) correspond to the end point of the hyperbolic 
region at cr = (k = 0), which is isolated in phase space {Qu,Iu}- Eq. (12) for the periodic orbits at 
this a can be solved analytically with respect to 0. Identifying the root 9{ri, nu/uy) with its definition 
(15) for hyperbolic orbits we realize that all short diameters M(l,2) bifurcate at the deformations, 

r?bif(M,n)^ \ = ^ (n-l,2,3,---,Af-l) (18) 

sin(7rn„/rit,) cos(n7r/2M) 

and at each bifurcation a new family of hyperbolic orbits M{nu, n^) with Mriy reflection points is 'born'. 
The second equation presents the same bifurcation points and shows explicitly that the bifurcation 
deformations 77bif are also identical to the corresponding divergences of the Gutzwiller amplitudes for 
short diameters, see Eq. (6.47) of Ref. [7]. Each of the emerging hyperbolic orbits Mi(M — n, 2M) with 
Ml repetitions and n from Eq. (18) coincides exactly with the corresponding short diameter MiM(l,2) 
repeated MiM times at the deformation 77bif. For instance, for the triply repeated short diameter 3(1,2) 
(Ml — 1,M = 3) there are two bifurcation points at the deformations rjbif ~ 2/\/3 and 2 where the 
primitive hyperbolic orbits (2,6) (n — 1) and (1,6) (n = 2), respectively, are born (see these orbits in 
Fig. 3.6 and discussion nearby in Ref. [19], also Ref. [14] and Fig. la there). However, the short diameters 
are isolated in the phase space of action-angle variables {Qu,Iu}- They emerge as terms of the periodic 
orbit sum which are additional to the families of hyperbolic tori (see a more detailed discussion below) . 
The contribution of the primitive short diameter 1(1,2) can be calculated by the original Gutzwiller trace 
formula, except near the circular shape [7, 19]. This formula will be improved near all bifurcation points 
(18) and the circular shape in §5.2. 

The long diameter orbits M(l,2) are also characterized by 2M reflection points and correspond to 
a specific isolated point in {Qmlu} space. They are related to the separatrix value cr = 1 (k = 1). 
Again, their amplitudes can be calculated with the standard Gutzwiller trace formula for isolated orbits, 
with the same exception near the symmetry-breaking point of the circular shape [7, 19] (see §5.3 for the 
improved solution in terms of Airy functions near this point). 

The limit of the circular disk (jj = 1) may in some sense also be considered as a (one-sided) bifurcation 
point: here the family of diameter orbits (with /C = 1) break into two isolated diameters with IC — and 
complicated hyperbolic orbit families (/C = 1) with n„ ^ oo, n„ ^ oo, and ri„ : n„ — > 1 : 2, when the 
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/„i 



(a) 




o=iif'^'> 






(b) 



Figure 2: Energy surface Iv{Iu) and curvature d^Iy/dl^ are drawn in the upper and lower panels, 
respectively, from Ref. [19]. 



deformation (ry > 1) is turned on. Inversely, the long and short diameters and hyperbolic orbits which 
have /C = and 1 in the ellipse, respectively, merge into the families of diameter orbits with IC — I 
as r] —f 1. The discontinuous change of /C at 77 = 1 is accompanied by a divergence of the diametric 
amplitudes in the standard SPM. This is the symmetry breaking problem discussed in the introduction 
and below in §5.2 and §5.3. 

Figure 2(a) shows the energy surface in action space, in the form of the curve ly — /t,(e,/„) at 
fixed energy s. Specific primitive orbits (with M — 1) are illustrated, with the arrows pointing to the 
corresponding stationary points /*: the short diameter (at /* = or cr = 0, with G* — 0, tt), the 
'butterfly' (or 'bow-tie') orbit, the long diameter (at /* = /„ , with cr = 1 and G* = ±7r/2), the 
rhomboidal orbits with 4 reflections, and the 'creeping' orbit (at /* = lu ) as the limit of a 'whispering- 
gallery' mode with number of reflections riy ~ 00 and winding number n„ = 1. The limits to the 
separatrix correspond to infinite values of Uy and n„ for hyperbolic or elliptic orbits with the ratio 
nu/uy going to 1/2 from either side (see also Ref. [14]). We use the same notation for both short and 
long diameters in terms of the integers n^, Uy and M like for the elliptic and hyperbolic one-parametric 
families, specifying them also by the stationary points in the phase space variables a (or /„) for all orbits 
and G„ for the isolated ones if necessary. 



3.4 Curvature 

A key quantity in the semiclassical theory in terms of the action-angle variables is the curvature K of 
the energy surface 

(tyj 



K 



dl?. 



da^ uj„ da^ 



Y — 
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The partial derivatives appearing on the right-hand side above are given in Appendix A. Figure 2(b) 
shows K versus /„. In the limit a ^ one finds the curvature for the twice repeated short diameters 
considered as primitive orbits [19]. For our definition of the (non-repeated) primitive orbits, one has the 
curvature Kg larger by a factor 2, i.e., 

which is finite and negative for all deformations. K stays negative for the entire hyperbolic part < a < 1 
of the curve, whereas it is positive for the elliptic part 1 < cr < o'er. At the critical points a = \ 
(separatrix) and at Ucr (creeping point), the curvature diverges. It tends to — oo as one approaches the 
separatrix from the hyperbolic side, and to -|-oo from the elliptic side. For a — > dcr it also tends to -l-oo. 



4 Phase Space Trace Formula in Action- Angle Variables 

4.1 Action-angle variables 

We now transform the phase space trace formula (3) from the usual phase space variables (r,p) to the 
angle-action variables (0,/). The latter are useful for integrable systems because the Hamiltonian H 
does not depend on the angle variables ©, i.e., H — H{I). For elliptic billiard one has from (3) 

<7sci(£) = ^^ ^^ E / ^Q« / '^Q" / < J < ^(^ - H{I'u,Q) 

X expj J [^„(J',/",i„) + (/" - /') • 0"] - ti^A , (21) 

where = {Qu,Qy} are the angles and I — {lu, Iv} the actions for the elliptic billiard defined in 
the previous section. For simplicity we omitted here and below the Jacobian pre-exponential factor of 
Eq. (3) because this Jacobian taken at the stationary points is always unity when we apply the improved 
stationary phase method for the calculation of the integral over phase space variables, as noted above. 

4.2 Stationary phase method and classical degeneracy 

As noted in the introduction, we emphasize that even for integrable systems the trace integral (21) 
is more general than the Poisson-sum trace formula which is the starting point of Refs. [4, 5] for the 
semiclassical derivations. These two trace formulae become identical when we assume that the phase of 
the exponent also does not depend on the angle variables ©, like the Hamiltonian. Then, the integral 
over angles in (21) simply gives (27r)" where n is the spatial dimension (n = 2 for the elliptic billiard), 
see Ref. [5]. In this case the stationary condition for all angle variables are identities in the 2tt interval. 
This is true for the contribution of the most degenerate classical orbits like elliptic and hyperbolic orbits 
with /C = 1 in the elliptic billiard. For the case of orbits with smaller degeneracy like the isolated 
diameters (/C — 0) in the elliptic billiard, the exponent phase is a strongly dependent function of some 
angles with definite discrete stationary points. We therefore need to integrate over such angles by the 
standard or improved SPM. Other examples are the equatorial orbits (JC — 1) and diameters along the 
symmetry axis (separatrix with /C = 0) in the spheroidal cavity {n — 3), the degeneracy parameters of 
which are smaller than the largest possible value /C = /Cmax — 2 for the elliptic and hyperbolic orbits in 
the meridian plane, or for 3-dimensional orbits. We have a similar situation also for the diameters with 
/C = 2 in the spherical cavity (/Cmax — 3), orbits along the symmetry axis for axially-symmetric cavities, 
and so on. Thus, the stationary conditions with respect to the angle variables for orbits with smaller 
degeneracies are not identities. Moreover, the stationary points in the cases mentioned above occupy 
subspaces of the phase space which are isolated in the rational tori that lead to separate contributions 
to the trace formula, except for the most degenerate orbit families, as we shall see below for the case of 
the elliptic billiard. 
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4.3 Stationary phase conditions 

Wc first take the integral over I^ in Eq. (21) exactly. Due to the energy conserving (5- function, we are 
left with the integrals over angles O", 6" and action I'^: 

9s.ie) = ^ReY^Jde:Jde:JdIi^^ 



X exp 



n 



(22) 



S^{I',I",t^)^- j ^ dl ■&{!). (23) 



We first write down the stationary phase equation for /, 

'ds^{i',r,t^) 



dlL 



-e'l = e'^-e: = 2^m„, (24) 



where M„ is an integer. The star means that we take the quantities at the stationary point /^ = lu- We 
now use the Legendre transformation (6), which reads 

SM',I",t^) + (I" - I') ■ 0" = S^i&",&\s) - I' ■ (©" - 0'), (25) 



S'„ (©',©",£)== / d&-I{&) 
J&' 



With the use of this transformation, the stationary phase conditions for angles 0„ and Q-^ are written 
as 

For the following derivations we have to decide which stationary phase conditions from Eqs. (24) 
and (26) are identities for the finite volume of the phase-space tori and which are equations for the 
isolated stationary points. For doing this, we have to calculate separately the contributions from the 
most degenerate (elliptic and hyperbolic) families [K, — 1) to the improved trace formula and those 
from diameters in the elliptic billiard. These two contributions are different with respect to the above 
mentioned decision concerning the integration over the angles 0. After the integration over one of the 
angle variables, say Qv, corresponding to the identity in the stationary phase conditions (26) due to 
an invariance of the action along the periodic orbit in Eq. (22), one gets Eq. (7) of Ref. [30] derived 
earlier by Bruno [38]. So, we get the result of Refs. [30, 38] within periodic orbit theory. Our phase- 
space trace formula (3) is more general because it can be applied for more exact calculations of the level 
density, without use of the stationary phase conditions like Eqs. (26), in terms of closed (periodic and 
non-periodic) orbits. 

Note that we have separate contributions coming from each kind of families and isolated orbits even 
near the bifurcation points (18) where we have the end point. Taking the deformation at a small distance 
from Tybif, we are left with two separate close stationary points and then use the Maslov-Fedoryuk 
theory [39, 40, 41, 42] like for caustic and turning points. Finally, after the integration by the improved 
stationary phase method, we look at the limit rj -^ r^bif to the bifurcation point. In particular, this idea 
of Maslov and Fedoryuk was applied in Appendix B for the calculation of the contribution of the long 
diameter at the separatrix. 

5 Trace Formulas for the Elliptic Billiard 

5.1 Elliptic and hyperbolic orbit families (/C = n — 1 = 1) 

Each family of elliptic or hyperbolic orbits with a common action occupies a two-dimensional finite area 
in the elliptic billiard. In this case, the stationary conditions (26) for the integration over the angle 
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variables 0„ and O^, become identities, since the integrand does not depend on the angle variables, and 
we have the conservation of the action variable I'^ = /" = /« fulfilled identically along each classical 
trajectory a. Taking the integrals over gives a factor (27r)^, and we are left with the Poisson-sum 
trace formula like in Refs. [4, 5]: 



5sci(£) 



1 



■Re 



M 



" M 



dl 6{e-H (I)) cxp 



27ri 



-MI- ivM 



din 



1 

\uj„ 



exp 



—-M ■ I - ivM 



(27) 



Here M = {MmMy) are integers which correspond to those in Eq. (14). Next we transform the inte- 
gration variable in the last expression of Eq. (27) from /„ to a defined by (9). Thus, the level density 
component SgscX.i related to the one-parameter families can be written as a sum of contributions from 
the hyperbolic ('5.9gJi\(e)) and the elliptic {Sg^^[-^{£)) parts of the tori. Their sum is 



5gsc\,i{e) 



1 



i:eopR? 



Re^ 



M 



1 

'^« Jo 



dalM -TT- exp 

0(J 



27ri 



M-/(cr) -iUM 



(28) 



where £q ~ fi /(2m_R^), /(cr) are the actions defined by Eqs. (11), Im arc the 'lengths' of the primitive 
orbits with ikf = 1 given by 

27rn^,p 



I 



M 



mujy 
2nyb sin( 



E(0,k) 



F(|,«) 



E(f,K) 



cot 



61^1 - k2 sin^ e 



(29) 



and 0(a) and k((t) are defined by Eq. (15). The 'lengths' become the true lengths of the corresponding 
periodic orbits when they are taken at a equal to the real positive roots of Eq. (12) inside the integration 
range. For other values of a, the 'lengths' are nothing else than the functions (29) introduced in place 
of Lu-u for convenience. The integration range from the bifurcation point cr = to the separatrix CTs = 1 
covers the contributions of all hyperbolic orbits. The remaining part of Eq. (28) from cr = 1 to the 
creeping value CTcr gives the contributions from the elliptic tori. 

As we shall see below, the choice of cr as the integration variable significantly improves the precision 
of the SPM. We hence apply the stationary condition (24) for the phase in the integrands of Eq. (28) 
with respect to a rather than to /„. With Eqs. (15), this condition becomes identical to Eq. (12) and 
determines the stationary phase point a' = a" = a* related to /^ = /" = /*. We used here the 
conservation of a (or the additional integral of motion 1+1-) along the periodic orbit. We now expand 
the phase up to second order. 



5„(J',/",t„) + (/" - /') • 0" = 27rM • / = Sfi{e) + \^j\{g - o*f, 
where Sfi is the action along the periodic orbit /3 determined by Eq. (12), 

Ssiie) = 2TrM{nuIu{<y*) + nyly{a*)), 
and ji is the Jacobian stability factor with respect to cr along the energy surface: 



4 = 






2ttM 



d^l, 






It is related to the curvature Kp (19) of the energy surface by 



J^ = 27: MuyKfj 



diu 
da 



27rMn„e \K^ 



/3| 



CT=<T*,/3 



dIu 
da 



(30) 



(31) 



(32) 



(33) 



cT=a*,f3 



where e = +1 for elliptic orbits and e = —1 for hyperbolic orbits. We substitute now the expansion 
(30) and take the pre-cxponential factor off the integral in Eq. (28). For the sake of simplicity, we only 
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consider the lowest order in the expansion ol the phase and the pre-exponential factor in Eq. (28) in 
the variable a, although higher-order expansions can in principle be used to improve the precision of the 
SPM. Thus, we are left with the integral from tr = to 1 for the hyperbolic orbits, and from ct = 1 to 
(Tcr for the elliptic orbits. 

When the stationary point a* is far from the limits of these intervals, one can extend the integration 
range from — oo to oo and get the result of the standard POT [4]. Near the bifurcation points (18) of the 
short diameter orbit (where the hyperbolic orbit families appear), however, the stationary point a* is 
close to zero. In this case we cannot extend the lower limit to — oo but have to take the integral exactly 
from (7 = 0. On the other hand, when the stationary point a* approaches the integration limits as (16) 
or (Tcr (17), hyperbolic or elliptic orbits with an increasing number of corners n^ appear. In these cases, 
too, we cannot extend the integration limits to ±oo. Taking the integral over a within the finite limits, 
we obtain a trace formula in terms of complex Fresnel functions or generalized error functions. The 
contributions of the one-parameter orbit families Sgsd.ii^) a-re then given in the form 

SgsciA^) ^'R.G^Af\e)exp ikLp ~ iiy^f^°^' . (34) 

Here, the sum is taken over both elliptic and hyperbolic orbit families, k — \/2me/h. The amplitude 
Aa of the orbit family /3 is given by 

^fl ^ = , '^ = erf (zl , , Z" „) ; (35) 

Lfj is the 'length' of the orbit family (29) corresponding to the stationary point a* [M ~ 1). We have 
introduced here the generalized error function erf(zi, Z2): 

erf(zi,Z2) = ^^ / dze^^ = erf(z2) — erf(zi), (36) 

erf(z) being the standard error function [47] with (complex) argument z. The complex quantities Zi ^ 
and Zi 2 ill (35) are given in terms of the Jacobian ji (32) and the stationary points a*: 





2^.1 = IHtt^ <i„ - '^* - 4,2 ^\-^[ <L -^*h (37) 



where a^^l^ and (Tmax are related to the integration limits by 

The phases J^i ° in (34) are related to the Maslov indices. They have a constant part 1^/3 which is 
independent of deformation rj and energy s. At deformations which are far enough from bifurcation 
points, such that the stationary points are far enough from the integration limits, we can determine 
this asymptotic part i^/s by transforming the error functions to Fresnel functions [47] with real limits 
and extending the integration limits to ±00. We hereby arrive at the amplitude Ag of the standard 
POT [4, 11, 44] 

^(1) = ^1^ (39) 

" eoTTkR^y/~eiM^nl\hKp\' 

and lyp is determined by the number of turning and caustic points as in the theory of Maslov and 
Fedoryuk [39, 40, 41, 42]. In terms of the numbers n^, and n„ and the repetition number M, it is given 

by 

3 

^^ = -2n^M for e = -|-1, 

TT 

i^iS^-{2nu + 2n^)M for e = -1. (40) 
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^From Eqs. (34), (35), and (40) we determine an extra contribution to the total phase ^"1° 



(tot) 



,^^;°'\v, kR) ^y!i-\e- arg {crl (4,i' ^.s) } (41) 



that analytically connects the asymptotic values Vf) and depends on the energy through kR. The final 
result (41) for the total phase depends also on the deformation parameter r\. 

Note that a* is negative for r\ <r7bif- In the derivation of Eqs. (34) and (35), we have changed the 

integration variable from cr to z = W— iej j!||/(2?i)(cr — a*^ in order to transfer the kR and r) dependence 
of the integrand to the limits of the complex generalized error functions (36) . Note also that our energy 
and deformation dependent phase Vq° is essentially different from Ref. [26] and much simpler in its 
analytical structure. Different from Refs. [26, 29], we did not use any assumption concerning a smoothness 
of the phase. Our solution is regular at the separatrix and creeping points, at all bifurcations points and 
in the circular disk limit. We easily get the correct circular disk limit [46] and the Berry- Tabor result [4] 
for larger deformations far from the bifurcations. 

Equations (34), (35) and (41) represent one of our central results concerning the contributions of the 
degenerate orbit families (/C = 1), that simultaneously solves the symmetry-breaking problem for both 
hyperbolic and elliptic orbits: near r; = 1 and other bifurcation points for all hyperbolic orbits, and near 
the separatrix Gs and the 'creeping' point o^t for all elliptic orbits. The additional contributions of the 
isolated orbits (/C = 0) will be derived in the two following subsections. 

Formally, our result (34) coincides with the first main term of the Berry- Tabor trace formula, see 
Eq. (24) of Ref. [4] , using the simplest way of the expansions near the stationary point instead of a more 
general and more complicated mapping procedure. The next two terms of their formula, being of higher 
order in v^, can be obtained by accounting for the linear term in the expansion of the pre-exponential 
factor over a — a* . They were neglected in our approach because we are interested here only in the main 
term of the SPM expansion, in order to get the simplest possible solution of the bifurcation problem. 
With the higher-order corrections, we should take into account that the ratio of the contribution of the 
linear term to the zero-order term of the amplitude is of the same order as the relative contribution of 
the next order (cubic) term in the expansion of the phase. For a consistent treatment of the level density 
in the semiclassical asymptotic approximation kR ^ 1, one would have to collect both corrections. 

5.2 Short diametric orbits (/C = 0) 

For the contribution of the isolated (/C = 0) diameters, only one of the two stationary phase conditions 
(26) corresponding to the O^, variable is an identity. The other one for 0„ is a nontrivial equation for 
the discrete number of the stationary points which differ by integer multiple of tt. Indeed, due to the 
integrability of motion in the elliptic billiard one has 

e„ = c^„t + el°) , e, = c^,f + e(,°) , (42) 

where 0^^) is the initial angle at t = 0. Since the frequency u;„ in Eqs. (42) is zero for short diameters, 
for instance, there is no room for an identity in the stationary phase condition for the variable 0„ in 
Eq. (22). Hence, the Poisson-sum trace formula cannot be applied to get the contribution from the short 
diameters unlike in the derivations in Ref. [24]. The stationary points for the integration in Eq. (22) 
over angle 0„ for the short diameters are constants 8* = ttM for M = 0, ±1, . . .. Due to the periodicity 
of the angle variable with the period 27r we really need to deal with the two stationary points 8* = 
and TT in the integration interval from — tt to tt over the angle 0„ in Eq. (22). We can then reduce the 
initial integration interval for angle variable 0„ to the region from — 7r/2 to 7r/2 taking into account 
the integration over other angles (related to the motion along the same periodic orbit in the opposite 
direction) by the factor 2 (due to the time reversal invariancc of the Hamiltonian) . Within this reduced 
integration interval, only one stationary point 0* = must be taken into account in the calculation by 
the improved stationary phase method. 

For the other variable O^ for the short diameters we have identity in the corresponding equation from 
Eq. (26). The integrand in (22) is independent of the variable 0t, and the integral gives simply 27r. Thus, 
the integrand for the contribution of the short diameters essentially depends only on 8„ and possesses 
the relevant stationary points. When we take this integral by the SSPM we get immediately Gutzwiller's 
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result for the short diameters with his stabiUty factor in the denominator. This stabihty factor is zero at 
the bifurcation points. We shall below get the short diameter term improved at the bifurcation points. 
For this purpose we shall first follow the same method in the integration over 0„ and /„ as we did in the 
integration over /„ for elliptic and hyperbolic orbits with highest degeneracies. The integration interval 
over /„ for the contribution of the short diameters is also finite from to the maximal "creeping" value 
lu (17) which corresponds to the region of the cr variable < cr < o'er- 

Thus, for the short diameters, we use the stationary condition for the angle variable 0„ and expand 
the phase of exponent in Eq. (22) about the short diameter, 

Sa = Ssm{£) + ■^Jsm'^u^ (43) 

with Ssm{£) being the action along the short diameter, Ssm{£) = 4:p{e)aM and 0* = 0. J^ is the 
Jacobian corresponding to the second variation of the action Sa with respect to the angle variable 0^, 

(44) 

according to Eq. (25). The Jacobian jjj^^ is expressed in terms of the diametric curvature Kg (20) and 
Gutzwiller's stability factor F^m, 

as'^ ^" ^^" =4sin2[Marccos(2rr^-l)], (45) 

which is independent of the choice of the phase space variables 

T-L IP T(e) FsM ,.„■. 



where 



7(e) 



dl' 



\^^uj sM 

and Ks is the short diametric curvature given by Eq. (20) (e = —1). In the second equality of Eq. (46) 
we used a simple relation between the Jacobians J^^ , JI and Kg- This relation follows directly from 
their definitions and simple properties of the Jacobians: 

T(e) .11 
1^=:-1. (48) 

\ da I 

After the exact integration over Q^ in Eq. (22) which gives 27r as explained above, we substitute the 
expansion (43) of the action Sa and take the amplitude factor at the stationary point 0* = 0. We take 
the integral over 0„ within the finite range from — 7r/2 to 7r/2 which can be reduced more to the integral 
from to 7r/2 with the factor 2 due to the spatial symmetry in addition to the time reversibility factor 
2 mentioned above. Integrating over /„ as in the previous subsection, one finally gets 



'^5sci,o == f'-^ X! -^sM exp[ifcLsM - Wsm\ ■ (49) 



M 

Here, Lsm is the length of the diameter orbit, Lsm = 4Ma, 



■^sM — T"rrD2 / I rn I ^^^ (^sM,l'^sM,2) erf (^^Mai-^gj^^j) ) (50) 



eo7rfci?2 yi^ 



sM\ 



ZsM,i and ZsM,2 are defined by 

(51) 



7" — n 7" 



\ 



i 7" 



2n 
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^t^^A/f„, ^ _l / ^ fJ c-n/T ^ ,, TT/ ^ ^ , 



2.V1 = \ — ^e; = 0, zi-M. = \^ii^e: = ;-j^ii^. (52) 



'sM 



2n " ' ■''"''' V 2?i " 2 V 2?i 

For any finite deformation and sufficiently large kR, Eq. (50) is much simplified by using asymptotics 
for the first error function and one gets 

The constant part VgM of the Maslov phases in Eq. (49) is obtained in the same way as in the previous 
subsection, 

VsM^iT^M-^. (54) 

For deformations far from the bifurcation points, the level density iJ^s^i ('^^) asymptotically reduces 
to the standard Gutzwiller formula for isolated short diameters, [1, 2, 7] 

•^^SoC^) - J^. E 7^ -^(fc^^M - ^sm). (55) 

The total Maslov phase v)^^j for the diameter orbits is 



(tot) 



''^^''O ■ TTSokR^ "^ ^ ^2Mih\K, 



V 



M 
1/4 



arg jerf (2[^^, z|^^^^ j | - arg {erf (Zi^.m, ^^^m) } 

~ t^sM - arg {erf (Z^^^ji^, Z^^^^) } (56) 

for large kR. 

Near the bifurcation points where Fsm ^ 0, one gets from Eq. (49) the finite limit, 

(^) - - ^^ Re y , ^ erf (Zj^, „ Zj^ 2) e.(feL.M-..M) 

, RgX^ 1 g«(feL.M-^aM-V4)^ (57) 

Note that the two last terms in Eq. (24) of Ref. [4] are smaller than the above contribution (57) at the 
bifurcation deformations ?7bif (18) by the factor ^/kR. Therefore, these two terms are the next order 
semiclassical corrections and can be neglected compared to the term (57) obtained above. Moreover, the 
ISPM solution (49) is not related to the "diametric" part of the Poisson-sum trace formula (28) with 
n„ = 1, n„ = 2 as follows from the derivations in Ref. [24] (ai = 2, a2 = A = 2 in the notations of 
Ref. [24] applied for the short diameters in the elliptic billiard, ai = 2n„) (see a more detailed discussion 
below). Thus, our derivation is essentially different from that suggested earlier in Ref. [24] (where the 
last two terms in Eq. (24) of Ref. [4] are retained without considering the contribution (57)). 

Taking the limit of Eq. (57) for r] ^ 1 we obtain the same contribution of the diameters in the circular 
disk [46] as found from the "diametric" part of the Poisson-sum trace formula, 

'552i(e) = j=^= E -Jfj sin(fcisM - VsM + 7r/4). (58) 

£oV 27rfci? ^ VM 

The value at this limit is larger by the factor Vfci? than the standard Gutzwiller result for isolated orbits 
like at any other bifurcation points. 

5.3 Long diameters and the separatrix 

As shown in §2, the curvature K goes to infinity being positive from the right side and negative from 
the left side near the separatrix (tr = 1) with the same modulus, see Eqs. (87), (88) and Fig. 2(b). The 
derivation for short diameters of the previous section with the expansion of the action exponent phase to 
second order terms cannot be applied in this case. However, we note that the behaviour of the curvature 
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near the separatrix in the action /„ (or a) variable is similar to that for the eigenvalues of the matrix of 
the second derivatives of the action in the usual coordinate space near the turning points. One can thus 
apply the Maslov and Fedoryuk idea for the calculation of the Maslov indices, see Refs. [39, 40, 41, 42]. 
Following this idea we expand first the phase of exponent in Eq. (21) with respect to the action /„ 
taking into account the next third order terms, see Eq. (89) in Appendix B. We then use the linear 
transformation (97) to the new variable z to get the standard exponent in the integral representation of 
the Airy functions. Within this method we take a small first derivative (small parameter ci) and the large 
second derivative (curvature) in the cubic polynomial expansions (89) taking a on the small distance from 
the separatrix a = 1. After some algebraic transformations wc obtain Eq. (100) in Appendix B for the 
limit a ^ 1. Note that a similar idea which we used here considering a near the singular separatrix point 
a = 1 and only finally, after the calculation of the integrals, put cr — > 1 was applied in the derivations of 
the separate contributions of the hyperbolic orbit family and short diameters to the periodic orbit sum, 
as mentioned above. 

For the angle integral in Eq. (100) we use the same Maslov-Fedoryuk method [39, 40, 41, 42] applied 
for the caustic case. As result, one obtains (see Appendix B) 



^ [Ai(-w7|| ) + I Gi{-w\\ )] 



--2''2 



\ 

Here, the complete and incomplete Airy (or Gairy) functions with one and three arguments (Eq. (102)) 
are used in line of the definitions in Refs. [47, 48], see also Appendix B for the definitions of all other 
quantities. 

For large kR^Jrf — 1, near the separatrix cr ^ 1 the parameter wj^ is negligible in Eq. (105) for the 
limits Z^ij^ and Z^if^i and the integration range can be extended from to oo. The incomplete Airy 
integrals in Eq. (59) approach the complete ones and the asymptotics of all Airy functions like Ai(— w) 
and Gi(— w) is now applied [47]. Finally, we get asymptotically the standard Gutzwiller result for the 
isolated diameters, [1, 2, 7] 



X [Ai(-W|| ) + i Gi(-W|| )] [Ai(-u;_L) + i Gi{-wi_)\ 

where Fim is the Gutzwiller stability factor for long diameters, 

FiM = -4sinh2 [M arccosh(2772 - 1)] , (61) 

i^iM = 37rM-|. (62) 

In the second equation we used the asymptotics of the Ai(— w) and Gi(— w) functions [47]. Wc found also 
the constant part vim of the phase by using the Maslov-Fedoryuk theory. The deformation and energy 
dependent Maslov phases are determined by the additional phases in the exponent and the argument of 
the product of the square brackets in (59) through complex combinations of the Airy and Gairy functions 
and their arguments. 

In the circular shape limit both the upper and the lower limits of the incomplete Airy functions in 
Eq. (59) tend to zero and the angle integral has the finite limit 7r/2 because Cj, c^ and w±^ vanish, 
see Appendix B. With this, the other factors near the separatrix ct — > 1 ensure that the amplitudes for 
long diameters diminish because W|| (99) vanishes at the separatrix, see also Ref. [47]. Namely, the long 
diameter contribution becomes zero in the circular shape limit. 
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Figure 3: Illustration of the caustic method for evaluating the stability factor Ja in Eq. (64) for the 
closed two-reflection orbit "co2". The deflection angle 66' at the initial point O(r'), variation Sy" of the 
final point 0'{r") with respect to O and the coordinate system {x,y) are shown: Thick solid lines and 
dashed lines show the hyperbolic orbit "co2" and the perturbed orbit, respectively; the thin solid curve 
indicates the orbit-length invariant hyperbola confocal to the boundary. 

Thus, for deformations far from the bifurcations, the results (49) and (59) of the ISPM reduce to 
the standard Gutzwillcr formula. In the circular disk limit the improved short diameter density (49) 
approaches continuously the diametric contribution to the circular disk density, while the long diameter 
(separatrix) contribution diminishes. Note that our ISPM solution (59) for the unstable long diameters 
is not related to the Poisson-sum trace formula (27), in particular, with its "diametric" part because 
of the existence of the isolated stationary points for the angle variable O^ like for the short diameters. 
Moreover, the uniform approximation Eq. (24) of Ref. [4] is singular at the separatrix because of the 
divergence of the curvature Ki for cr ^ 1, as noted in Ref. [26]. However, instead of the suggestion 
of Ref. [26] to use the continuation of the WKB approach to the complex plane we applied simpler 
Maslov-Fedoryuk method [39, 40, 41, 42] and got the analytical dependence of the Maslov phase on the 
deformation and energy through the exponent phase and complex arguments of the Airy functions as 
well as their complex summations. 

5.4 Closed orbits and the circular disk limit 

To get a more exact solution for the diameter contribution to the level density and check the precision 
of the ISPM, we come back to the initial trace formula Eq. (2) before application of the ISPM for 
the calculation of this trace. ^ For this purpose we shall take exactly the trace integral (2) in suitable 
variables. This is the trace formula in terms of the sum over all closed (periodic and non-periodic) orbits 
a, 

dxdy ■ /, T \ fnn\ 

sin(kLa — i^a), (63) 



<55sc>(.) = 2(2.n)-/^^5:/^ 



a J \/Ja{x,y) 

where Ja{x,y) is the stability factor defined through the Jacobian i/a (p'^q , f "e) by 

J'„(p'i„,r"£) = — — ^ =— — ^. (64) 



^Eq. (1) can be obtained also from the phase space trace formula Eq. (3) taking the integral over two components of 
the momentum p' along the energy surface by the stationary phase method. 
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Figure 4: Closed non-periodic two-reflection orbits with the elliptic and hyperbolic caustics at the initial 
point 0{x, y) are shown by thin and thick solid lines, respectively, for the deformation 77 = 1.05 (left-hand 
side) and 1.2 (right-hand side): O is the common vertex of both triangular orbits; dashed curves indicate 
the orbit-length invariant ellipse and hyperbola crossing the initial point. The hyperbolic orbit is close 
to the diameter of the circular shape for small deformations. 



Here the deflection 5y" of the final path point in the perpendicular direction of the local Cartesian system 
{x, y) comes from the angle variation 56' of the initial momentum, [11, 46] see Fig. 3. 

We shall then simplify the trace formula (63) taking the contribution of the main shortest closed 
orbits a with the two reflection points denoted below by the index "co2" as an example. For arbitrary 
point (x, y) inside the elliptic billiard one can find such orbits "co2" which are the triangles with the two 
vertices at the elliptic boundary and one vertex at the point {x,y), see Fig. 4. There are two kind of 
such orbits. For any point (x, y) we can plot the hyperbola and ellipse confocal to the boundary, which 
are the orbit-length invariant curves. Indeed, moving the initial point {x, y) along such hyperbola (or 
ellipse) we have the one-parametric family of the triangle- like orbits with the same action (/C = 1). We 
shall call them for short the hyperbolic and elliptic "co2" orbits, respectively. 

For the calculation of the trace integral (63) it is convenient to use the elliptic coordinates (u,w), 
(7). After this coordinate transformation, we can take the sine function of the action off the v or u 
integration for the hyperbolic or elliptic "co2" orbits, respectively, because of independence of the action 
on the corresponding elliptic coordinate. Finally, one obtains from Eq. (63) 



&ie) 



2{2TTny 



■^li°rnC^ f dusm{kL\ico2iu) ^ Vhco2) dv(smh v + cos^ u) 



Vp 



^/Jhco2{x{u,v),y{u,v)) 



(65) 



for the contribution from the hyperbolic "co2" orbits (hco2), and a similar equation for the elliptic "co2" 
orbits. An explicit expression for the stabihty factor Jco2{x,y) evaluated by the caustic method [11] is 
presented in Appendix C. 

Note that the hyperbolic "co2" orbits with the initial point (x, y) reduce to the disk diameters crossing 
the same point in the circular disk limit, sec Fig. 4. The stability factor Jhco2{x,y), (108), turns into the 
analytical circular disk expression of Rcf. [46]. The circular disk limit of the level density (65) coincides 
with the diameter contribution '^.9sci i(^)' (5^)' ^^ shown in Fig. 5(a). The opposite limit of (65) far from 
the bifurcations is the Gutzwiller SPM for the short and long isolated diameters, see Fig. 5(b). The 
contribution of the elliptic "co2" is negligibly small everywhere, and vanishes at the circular disk shape 
as next order h corrections. 
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Figure 5: (a) Convergence to the circular shape hniit: The contribution of the closed two-reflection orbits 
of the hyperbolic type "hco2" (see Fig. 4) to the level density Sg{kR) is shown by a solid line for the 
deformation i] — 1.005, while Gutzwiller's trace formula (SSPM) for isolated diameters and circular disk 
trace formula are indicated by dotted and dashed lines, respectively. The dashed line overlaps with the 
solid line, so that it cannot be distinguished from the latter, (b) Convergence to the Gutzwiller trace 
formula for rj = 1.1. Notations are the same as in (a). 

6 Level Density, Shell Energy and Averaging 

6.1 Total level density 

The total semiclassical POT density can be written as the sum over all periodic orbit families considered 
in the previous section. 



Sgsciis) = Sg^^^A^) + Sg%{e) + Sgil,{e) = ^ Sgi^J{e), 

13 



(66) 



where the first term is the contribution (34) from the elliptic and hyperbolic orbits. The second and 
third terms are the contributions from the short (49) and the long (59) diameters, respectively. Near 
the circular limit, the last two terms for one period (Af = 1) can be replaced by the contribution of the 
hyperbolic "co2" orbits (65) to get a more precise semiclassical result. 



6.2 Semiclassical shell energy 



,(/3)^ 



The shell-correction energy 6E can be expressed in terms of the oscillating part (^gg^fj (e) of the semiclas- 
sical level density as [6, 7, 11] 



--i:(^)^*.2,'M. «-/■ 



deg{e). 



Here, tfj is the time of the motion along the periodic orbit [3 (including its repetitions), 

27rM/3 



tp = MpTp 



nf3 



(67) 



(68) 



where T^ is the period of primitive orbit with the Fermi energy ep, M^ the repetition number, f2^ the 
frequency, and N the particle number. Note that we have taken into account the spin degeneracy factor 
2 in (67). 

The semiclassical representation of shell-correction energy (67) differs from that of Sg only by a factor 
{h/tpY' = {h kp/mLf^)'^, which suppresses contributions from longer orbits. Thus short periodic orbits 
play dominant roles in determining the shell-correction energy. 
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6.3 Average level density 

For the purpose of the presentation of the level density improved at the bifurcation points we need to 
consider a level density averaged slightly, thus avoiding the convergence problems that usually arise when 
one is interested in a full semiclassical quantization. 

The averaging is done by folding the level density with a Gaussian of width T: 

g^{e) = -^j^ de'g{e')e-^^) . (69) 

The choice of the Gaussian form of the averaging function is immaterial and guided only by mathematical 
simplicity. For cavities it is convenient to use also the level density defined as a function of kR averaged 
with a Gaussian of width 7: 



gJkR) = ^^ d{k' R) g{k' R) e \ -< I , (70) 

y/T^l J -00 



where 



g{kR) = V 5{{k - ki)R) = 2kReo V S{e - e,) = 2kReog{e), (71) 



So = ^ /2mR^ and the dimensionless parameter 7 is related to F by 

F - 27^^^. (72) 

Applying the averaging procedure defined above to the semiclassical level density (66), one gets [3, 
46,11] 



S9r.Ms)-E^9i'Jie)e-("^') = E ^5L1^(=) ^"^"" ^ 



(73) 

The latter equation is written specifically for billiard problems in terms of the orbit length Lp (in units 
of a typical length scale R) and 7. The averaging yields an exponential decrease of the amplitudes with 
increasing Lp and/or 7. As shown in Ref. [11], for 7 of the order of unity, all longer paths are strongly 
damped and only the shortest periodic orbits contribute to the oscillating part of the level density, 
yielding its gross-shell structure. For a study of the bifurcation phenomenon, however, we need smaller 
values of 7. 

Finally, we should note that the higher the degeneracy of an orbit, the larger the volume occupied by 
the orbit family in the phase space and also, the shorter its length, the more important its contribution 
to the average level density. 

7 Quantum Elliptic Billiard 

7.1 Numerical method for the spectrum calculation 

Single-particle energies Si of a particle of mass m moving freely inside the elliptic boundary v < Vh can be 
obtained by a number of numerical methods. Following the procedure employed in previous works [18, 20] 
by some of the present authors, one can expand the deformed single-particle wave functions \E'(r, 9) into 
a circular basis with well-defined orbital angular momentum I: 

(e) (o) 

*|++^(r,0)=EAiJi(fc,r)cos(/0), *(-+)(r,0) =EB,JKMsin(Z0), 
i=o 1=1 

^l+-\r,e)^J2'^iMhr)cosiW), ^\—\r,e)=J2BiMhr)sm{W), (74) 

1=1 1=2 



where Ji{x) are the cylindrical Bessel functions of the first kind, ki = ^/Qmel/ti, the superscripts (+-I-) 
etc. stand for parities with respect to reflections about the x and y axes, and the superscripts (e) and 
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Figure 6: Single-particle spectra (in unit of £o) for the elliptic billiard plotted as functions of the defor- 
mation parameter 77. 



(o) indicate the sums with respect to even and odd /, respectively. The expansion coefficients Ai and Bi 
can be determined by applying Dirichlet boundary conditions. 

In the present analysis, we also employed, in addition to the above circular-wave decomposition 
method, the numerical procedure based on a rather standard approach, the separation of the Schrodinger 
equation in the elliptic coordinate system [26, 52, 53]. In terms of the elliptic coordinates (7), the 
Schrodinger equation can be written as 



d 



d 



^^^a? H^a? 



d 



d 



Vi^^TH vT^^TTi i^{L4>) 



+ 



2me,e{e 



^(e,0) = o. 



(75) 



where ^ = coshv and (p — cosu. Following Ref. [52], this equation can be separated into two ordinary 
differential equations by assuming ^(^, (/)) = i?(^)S'(0). The functions R and S are solutions of the 
ordinary differential equations 



(1- 



l)^^^#^+e^^%^-[A.-c^^^]i?.(c,0 



de 



2.(PSl{c,( 



r- 



d^ 

dSi(c,(j)) 



0, 



+ [\i-c'4>"\Si{c,4>)^Q, 



(76) 



where A/ is the separation constant and c = C,\/2mei/fi for ^ < ^6 = coshw^. The internal radial functions 
Ri{c,S,) are expanded in terms of Bessel functions of the first kind. The expansion coefficients and the 
separation constant A; can be determined from the three-term recurrence relations found in various 
references [47, 52, 53, 54]. 

By imposing usual boundary conditions on the radial wave functions, i.e., Ri{c,S,b) = 0, one finds 
the eigenenergies Si. All eigenvalues up to kR sa 40 with the coordinate-transformation method can be 
calculated numerically in matter of minutes without overlooking solutions near level crossings, and hence 
the procedure is certainly effective for the present model. The results obtained from both numerical 
procedures were carefully compared and found to achieve a nice convergence. 

In Fig. 6 the deformation dependence of the single-particle energies for the elliptic billiard is presented. 
At the circular limit, the familiar shell gaps are clearly observed, while different shell gaps start to 
develop at higher deformations. Below we identify the semiclassical origin of these shell structures at 
higher deformations. 



7.2 Strutinsky's smoothed level densities and shell energies 

With the aid of the Strutinsky averaging procedure, [57] clear oscillatory patterns of the coarse-grained 
level density emerge as shown in Fig. 7, where (a) and (b) are obtained with the Gaussian smoothing 
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Figure 7: Coarse-grained level densities with the Gaussian smoothing parameter 7 = 0.3 (a) and 0.64 
(b). 
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Figure 8: Shell structure energy SE (in unit of £o) plotted as a function of both deformation 77 and 
particle number N. 

parameter 7 (defined by (72)) of 0.30 and 0.64, respectively. As clearly seen from these figures, the 
choice of a Gaussian smoothing parameter 7 is rather crucial for properly identifying the coarse-grained 
level density, and hence the contribution of classical periodic orbits. At the circular limit rj — 1.0, both 
Gaussian-smoothed level densities show similar oscillations, whereas the shell gaps for 7 = 0.64 start to 
collapse with increasing deformation. In particular for deformations rj larger than 1.5 the strong shell 
patterns cease to exist for the case of 7 = 0.64, while for 7 = 0.3 the appreciable effects still remain and 
show more oscillations as deformation increases. 

In the semiclassical picture, for a given value of 7 the contributions from only those periodic orbits 
of length up to £max ~ ttR/j can be considered. In this context, it is important to locate the actual 
shell-energy minima, irrespective of the choice of a Gaussian smoothing parameter. 

In terms of the particle number TV one can also obtain the shell-correction energy 6E defined as 
the difference between the sum of single-particle energies of N lowest levels (taking the spin-degeneracy 
factor 2 into account) and the Strutinsky averaged energies, i.e.. 



N 



£f 



5E^^e^-E, E = 2 de' e' g{e'), (77) 

with the Fermi energy ip satisfying 

7V = 2 / "^ de'g{e'). (78) 



Figure 8 illustrates the oscillating pattern of the shell-correction energy SE as function of both 
deformation rj and particle number N. It is clear from the figure that the distance between major 
shell gaps shrink with increasing deformation. In the considered range of deformation it is found that 
the actual magic numbers determined by the above procedure cannot be reproduced with the choice of 
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Figure 9: Smoothed shell-correction energies for 77 = 1.5 with Gaussian smoothing parameter 7 — 0.3 
(dashed line) and 0.6 (dotted line). Those without the smoothing are plotted by a solid line. 

7 — 0.64, whereas the value of 7 = 0.3 is small enough not to demolish but still large enough to keep the 
actual coarse-grained shell structure. It is explicitly shown in Fig. 9, where the shell-correction energies 
are now calculated by applying a Gaussian smoothing parameters of 7 = 0.3 and 0.64, for the case of 
77 — 1.5 as an example. In this case, the actual magic numbers are found to be ... , 16, 22, 30, 38, 52, . . ., 
which exactly coincide with those of 7 = 0.3, while those calculated with 7 = 0.64 show larger oscillations 
missing . . . , 16, 30, . . .. The same is true for other deformations considered in this paper. So, the coarse- 
grained shell structure obtained with 7 = 0.64 is too rough and therefore we adapt 7 ~ 0.3 to improve 
the precision of its description. 



7.3 Shell Structure and Fourier Spectra 

Equations of single-particle motion in the billiard are invariant with respect to the scaling transformation 



(r,p, t) -^ (r, ap, a t). The action integral Sp for a periodic orbit (3 is proportional to its length L^ 



Sfj{E = p^ /2m) = f dr ■ p = pLp = TikLp. 
Jp 



and the semiclassical trace formula for the level density is written as 

TT 



9sci{s) = g{e) + ^ Af3{kR) cos (kLp - -^p 



(79) 



(80) 



13 



where ^(e) denotes the smooth part corresponding to the contribution of zero-length orbits, Ap =^ \Ap\, 
fip the Maslov phase (the deformation and energy dependent phase of Eqs. (41) and (56) in our improved 
semiclassical approximation). As previously discussed, the stationary phase approximation employed in 
deriving the Gutzwiller trace formula breaks down at bifurcation points for stable periodic orbits, and 
consequently results in the divergence of the amplitudes Ap{kR) in Eq. (80), whereas in the present ISP 
treatment those amplitudes are smooth functions of both deformation and energy. 

In order to examine the classical-quantum correspondence on shell structure, one can perform the 
Fourier transform F{L) of the quantum level density g{e) with respect to the wave number k 



F{L) 



dk e 



-ikL 



9{sy 



2 Uc i 



2eoi?2 ^^ fc. 



-ikiL^ 2(^7) 



(81) 



which may be regarded as 'length spectrum' exhibiting peaks at lengths of individual periodic orbits. 
Here the Gaussian factor is imposed to smoothly cutoff the spectrum in the high-energy region. In 
numerical calculations, we use kc = fcmax/v2, fcmax being the maximum wave number included. The 
above method of taking Fourier transform of the quantum level density is known to be a powerful tool to 
investigate the role of classical periodic orbits in the appearance of shell fluctuations in quantum systems, 
and from such observations one can also extract the semiclassical contributions of individual periodic 
orbits. 
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Figure 10: Fourier transforms of the single-particle level density for the elliptic billiards with 77 = 1.0 
(a), 1.2 (b), 1.5 (c) and 1.7 (d). Some periodic orbits that correspond to peaks are illustrated. 



Fourier spectra for deformations ry = 1.0, 1.2, 1.5 and 1.7 are presented in Figs. 10(a), (b), (c) and (d), 
respectively. At axis ratio ry = 1.0, the diameter and elliptic orbits are found to be equally important. 
The fact that only those shorter periodic orbits mainly contribute to the gross-shell structure implies the 
significance of three classical periodic orbits at the circular limit, namely the diameter, triangular, and 
square shape orbits. As deformation increases, the Fourier amplitudes for triangular and rhombic orbits 
still exhibit fairly strong effects, while those for diameter orbits start to decline quickly and significant 
rearrangement can be observed. Especially at deformations 77 = 1.5 and 1.7, one can conclude, in 
addition to triangular and rhombic shape orbits, the gross-shell fluctuations are also governed by the 
(1,4) hyperbolic orbits bifurcated from the 2(1,2) short diameter orbit at the critical deformation rj — \/2- 

Figure 11(a) demonstrates the deformation dependence of Fourier amplitudes calculated from the 
quantum single-particle spectra. Here the enhancement of peaks indicates a larger contribution from 
the corresponding classical periodic orbits (3 of length Lp to the shell structure. At the circular limit 
the system possesses the highest symmetry, and the breaking of this symmetry due to a small deviation 
of its shape results in the orbital bifurcation. With increasing deformation, the short diameter orbits 
with M repetitions M(l,2) also bifurcate and create hyperbolic orbits at the critical deformations r^bif 
given by Eq. (18). The length of those classical periodic orbits as a function of deformation can be 
calculated [14] as shown in Fig. 11(b). It is clearly seen from both figures that the bifurcations of stable 
periodic orbits give rise to an increase in the Fourier amplitudes. The remarkable enhancements seen in 
the figure exactly coincide with the corresponding lengths of the newly created hyperbolic orbits, and 
hence stress the importance of the orbital bifurcations. 

In this context, similar enhancements for the case of a spheroidal cavity at superdeformed shape were 
also reported in Ref. [21], where superdeformed shell structure is associated with bifurcations of periodic 
orbits with two repetitions on the equatorial plane. In the present work, particular attention is paid 
to investigate such correlations between bifurcations of stable periodic orbits and quantum level-density 
oscillations. 

In Fig. 12, Fourier peak heights for some of the important hyperbolic orbits, namely those bifurcated 
from the short diameter orbits of 2, 3, and 4 repetitions, 2(1,2), 3(1,2) and 4(1,2), are displayed as a 
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Figure 11: (a) Modulus of the Fourier amplitudes plotted as functions of both orbit-length L and deformation rj. (b) Lengths L of classical periodic orbits ^ 

calculated as function of deformation 77. Solid, dashed and dash-dotted lines are used for hyperbolic, short-diameter, and elliptic orbits, respectively. 
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Figure 12: Deformation dependence of Fourier peak heights for hyperbolic and short diametric orbits 
2(1,2), (1,4), 3(1,2), (1,6), 4(1,2) and (1,8). Sohd hncs are used for multiple traversals along the short 
diameter, M(l,2) with M = 1, 2, 3, while long-dashed, short-dashed and dotted lines are used for hyper- 
bolic (1,4), (1,6), (1,8) orbits, respectively. Open circles indicate the bifurcation points. 

function of deformation rj. Interestingly, the Fourier peaks for these newly created orbits show a rather 
universal deformation dependence, that is, their heights reach the maxima shortly after their bifurcation 
points and quickly decrease with increasing deformation. Such remarkable features were already seen in 
Fig. 8, where the shell valleys for rj approximately larger than 1.5 can be understood to vary along the 
constant-action lines S{k,r]) — const, of the (1,4) hyperbolic orbits, as explained below. 

Suppose some classical periodic orbits P of length L^ are the dominant components in the semiclassical 
trace formula for the oscillating level density, then the shell valley maxima/minima follows the constant- 
action lines Sj3{k,ri) = const, of those dominating classical periodic trajectories. Referring to Eq. (80), 
such lines are determined by 



kLi3~ -lJ.[j = (2n+l)7r, 



0,1,2,. 



(82) 



We shall now demonstrate the above dependence in Fig. 13(a), where the smoothed level densities 
are plotted on the k-rj plane. As indicated in Fig. 13(b), it is interesting to note that the shell valley 
structures seen in Fig. 13(a) can be described by the constant-action lines of three major periodic orbits; 
near the circular limit the shell valleys vary along those of elliptic (mainly triangular and rhombic) orbits; 
in the right-half region of Fig. 13(a) the influence of newly created (1,4) hyperbolic orbits is visible; the 
contribution of short diameter orbits are less pronounced but certainly non-negligible throughout the 
considered range of deformation. The equality, Eq. (82), indicates the inverse proportionality relation 
between the orbital length Lp and wave number k. As the length of a trajectory (3 increases, the values of 
k decrease and consequently the smoothed level densities show more oscillations. In particular, since the 
length of the (1,4) hyperbolic orbits gradually increases within rj k, v2— 1-7 and then slowly decreases for 
•q > 1.7, the corresponding constant-action lines also behave in the same manner. Such a tendency was 
already observed in Fig. 8, where the contribution from the (1,4) hyperbolic orbits to the shell energy 
6E is apparent in the region rj > 1.5, indicating the essential role of the orbital bifurcations in quantal 
shell formations. 



8 Comparison between Quantum and Semiclassical Calculations 

Figures 14,15,16 show the modulus of the complex amplitude for a few short orbits. The semiclassical 
amplitudes for the hyperbolic "butterfly" M(ri„,n^) = (1,4) and elliptic triangular (1,3) orbit families 
calculated by the ISPM are in good agreement with the exact calculation of the Poisson-sum trace integral 
(22), see Figs. 14 and 15, respectively. All ISPM amplitudes are continuous function of the deformation 
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Figure 13: (a) Smoothed level density plotted on the k-Tj plane, (b) Constant action lines on the k-rj plane for the elliptic (1, 3) orbit (dash-dotted lines), 
the primitive short diameter 1(1,2) orbit (short-dashed lines) and the hyperbolic (1,4) orbit (solid lines). 
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Figure 14: (a) Amplitude modulus A for bifurcating short diameter 2(1,2) and butterfly 1(1,4) orbits 
obtained by ISPM are shown by solid lines as functions of deformation parameter 77; standard results 
of the extended Gutzwiller periodic orbit theory (SSPM) are shown by short-dashed lines, (b) ISPM 
amplitudes for the butterfly orbit (solid line) are compared with exact calculation of the Poisson-sum 
trace formula (27) (dashed line marked by "Poisson") and SSPM of Berry and Tabor [4] (dotted line). 
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Figure 15: (a) Same as in Fig. 14(a) but for primitive short diameter 1(1, 2) and triangle 1(1, 3) orbits for 
smaller deformations, (b) Comparison of the amplitudes for 1(1,3) with exact calculations and SSPM. 
Notations are the same as in Fig. 14(b). 
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Figure 16: ISPM amplitude modulus (solid line) for the sum of short and long diameter 2(1,2) orbits 
is compared with the (n„ = 1, n« = 2, M = 2) part of the Poisson-sum trace formula (27) (long-dashed 
line) and the Gutzwiller SSPM (dotted line). 



through the bifurcation point r/ = v2- A remarkable enhancement of the butterfly amplitude is seen at 
the deformation 77 = 1.5-1.6 slightly on the right of the bifurcation point (see Fig. 14). 

The ISPM amplitude for the primitive short diameter 1(1,2) quickly approaches the Gutzwiller SSPM 
result as one goes away from the circular limit and, for larger deformations, its magnitude is relatively 
small compared with those of other orbits mentioned above (see Fig. 15). 

In Fig. 16 we compare the ISPM result with the modulus of "diametric" part of the Poisson-sum trace 
formula corresponding to riu — 1, n-^ — 2 and M — 2, which is regarded in Ref. [24] as to represent short 
and long diameters, as well as the standard Gutzwiller results. The ISPM amplitude for the bifurcating 
short diameter 2(1,2) has the two maxima; at the bifurcation deformation \/2, which is significantly 
larger than the butterfly and triangular amplitudes, and at the circular shape, see also Figs. 14 and 15. 
(Similar maxima at the circular shape appear for any short diameter orbit. The maximum for the short 
diameter 1(1,2) is the largest one, in particular, larger than for the triangular orbit, see Fig. 15(a).) As 
seen from Fig. 16, there is the same circular shape limit for the ISPM approach and the "diametric" 
part of the Poisson-sum trace formula, which is identical to the diameter family amplitude in the circular 
disk. 

Apparently, the behaviour of the ISPM amplitude for two repetitions of the short diameter 2(1,2) is 
essentially different from that of the "diametric" part of the Poisson-sum trace integral which exhibits no 
enhancement near the bifurcation point. Thus, the Poisson-sum trace formula (27) describes the families 
with maximum degeneracy like hyperbolic and elliptic orbits, rather than the isolated diameters. For the 
isolated orbits with smaller degeneracy like diameters in elliptic billiard the Poisson-sum trace formula 
cannot be applied because of the isolated stationary points for the angle 8„ variable. This is the reason 
for the agreement of the ISPM and SSPM asymptotics unlike for the "diametric" term of the Poisson-sum 
trace integral in Eq.(27). It means that the diameters cannot be included in the usual EBK rational 
torus quantization. However, the diameters could be included in a more general quantization rule in 
terms of the averaged ISPM level densities (66) in a similar way as pointed out in Refs. [9, 12]. 

We note a significant improvement of the ISPM results compared to the SSPM for a close to the 
separatrix value 1 and the creeping value CTcr (17). These cases might seem to be important only in the 
limit rj —f 00 when CTcr tends to unity. However, even for < rj ^2 we meet the situations where the 
stationary points are close to the critical points cr = 1 and a = o'er, so that we have to integrate within 
the finite limits. 

We compare in Fig. 17 the semiclassical level densities Sgsc\{kR) calculated by the ISPM with the 
quantum results for the averaging parameter 7 ~ 0.3. The results obtained by the ISPM are in good 
agreement with quantum results even near the bifurcation point v2, where the SSPM gives the divergent 
result due to the zeros of the stability factor Fsm for the short diameters 2(1,2). For the deformations 
like 1.2 and 1.7 far from the bifurcation, one obtains a fair agreement between the ISPM and the SSPM. 
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Figure 17: Quantum and semiclassical (ISPM) oscillating level densities 5g{kR) versus kR for several 
deformations. Averaging parameter 7 = 0.3, parameter of the Strutinsky's shell correction method 
7 = 2.0 and correction polynomial degree 2M. — 6 are used. 
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Figure 18: Oscillating level density 5g(kR) versus kR (left-hand side) and shell energy SE in units eo 
versus N^^"^ (right-hand side) for a small deformation 1.01. Solid and dotted lines indicate results of 
quantum and ISPM calculations, respectively. Parameters of the Strutinsky's shell correction method 
are the same as in Fig. 17. 



Figure 19: Quantum and ISPM shell energy SE (in unit of eo) are plotted by solid and dotted lines, 
respectively, as functions of N^^^. 
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Figure 18 shows a nice convergence of the ISPM results to those of the circular disk trace formula 
for 77 — > 1. This convergence is seen for any small deformation when the semiclassical parameter kR 
becomes sufficiently large. With the inclusion of the closed (periodic and non-periodic) hyperbolic 
orbit contribution, one gets even better agreement with the quantum densities near the circular disk 
shape. For deformations far from the circular shape (?7^1.1) and far from other bifurcation points, 
the contribution of the hyperbolic "co2" orbits approaches Gutzwiller's SSPM result for the isolated 
diameters, see Fig. 5(b). 

For the averaging parameter 7 — 0.64, we have good convergence of POT sums for the ISPM and 
SSPM with a few short periodic orbits with M < 1, n„ = 1 and Uy < 10. This is due to the damping 
factor in Eq. (73) which ensures the convergence of the POT sum. For smaller 7 ~ 0.3 we need more 
orbits with M < 2, n^ < 2 and n„ < 10. Note that for 7 = 0.3 we have much better agreement of the 
ISPM results with quantum mechanical calculations than in the case of SSPM for the deformations near 
the bifurcations including the transition to the circular shape, see Fig. 7. 

Figures 19 and 20 show nice agreements of the ISPM results for the shell-correction energies with the 
corresponding quantum results. Note that we substitute the exact Fermi energy ep into the semiclassical 
shell energy SE (67) by using the second equation for the particle number there and quantum level 
density like in Ref. [11]. This is important to get the correct behaviour of the shell-correction energy 
as a function of particle number N, as explained in Ref. [11]. It is evident from Fig. 20 that the nice 
agreement between the ISPM and quantum results in the strongly deformed region of ry > \/2 cannot be 
attained without including the contributions from bifurcating 2(1,2) and (1,4) orbits. 

In all our calculations wc used the semiclassical approximation improved at the bifurcation points 
which becomes better with increasing kR for all deformations including the bifurcation points. 

9 Conclusion 

The most essential new result of this paper in comparison to the Berry- Tabor theory are two additional 
terms (second and third ones in Eq. (66)) in the improved trace formula for the elliptic billiard. These 
two terms represent the contributions from the short and long diameters which are continuous functions 
through all bifurcation points. For deformations far from the bifurcation points, we obtain asymptotically 
the standard Gutzwiller result for the isolated diameters, and the correct trace formula for the diameters 
in spherical limit of the circular billiard. Our results for the hyperbolic and elliptic orbits improved near 
the bifurcation points are simpler than those suggested within the uniform approximation [4, 26] . 

With the use of our improved trace formula, we have demonstrated the importance of bifurcations of 
the repeated short diameter orbit for the emergence of shell structure at large deformations. 
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Figure 20: Shell energy maps SE drawn as function of N^^^ and deformation 77. (Left-hand side) Quantum results like in Figs. 16, 18 and 19: (Middle- 
and right-hand sides) Semiclassical ISPM results with and without taking into account the bifurcating orbits, respectively (see text). 
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Appendices 
A Curvatures 

The actions /„ and /„ given by Eq. (11) are expressed explicitly in terms of the elliptic integrals [47, 49]. 
For elliptic orbits one has 

TT 1 



lu = -C\/2TOecr E -, —j= 



ly = —(V2mea 
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For hyperbolic orbits, 
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Eqs. (83) and (84) may be regarded as equations for the energy surface e{Iu,Iv) written in terms of the 
parameter a for its elliptic and hyperbolic parts, respectively. 

The curvature K of the energy curve are obtained by differentiating Eqs. (83) and (84) with respect 
to the parameter a. In this way one gets Eq. (19) with the following derivatives for elliptic orbits, 
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For hyperbolic orbits. 
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With Eq. (85) we obtain the curvature Kp (19) for elliptic orbits as 
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B Separatrix 

Like for the ease of turning points [39, 40, 41, 42] one writes 



= ^0 +n^+ 3^ 



Here, 
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where the star means /^ — /" — /* . The asymptotic behaviour of the constants c,[ near the separatrix 
a K 1 was found from 
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The second equality in Eq. (89) was obtained by a hnear transformation with some constants a and /3, 

x = az + l3, a=(34)-^/^ /? = -4/(34), (97) 

41 - (co - ciC2/(3c3) + 2c3/(27c2))ll, t| = a[ci - c2/(3c3)]ll. (98) 

Near the stationary point for cr ^ 1, one has cl ^ and rj -^ — wn with the positive quantity 
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Using expansion (89) in Eq. (21) and taking the integral over angle 0" exactly, i.e. writing 2tt instead 
of this integral, one gets 
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Ai(— ui, zi, Z2) and Gi(— w, zi, Z2) are incomplete Airy and Gairy functions, [48] 
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(102) 



and Ai(— w) and Gi(— w) are the corresponding standard complete functions [47]. We used in the second 
equation of Eq. (100) that for any finite deformation rj and large kR near the separatrix {a -^ 1) one 
gets (see Eq. (99)) 
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Using an analogous expansion of the action tq in Eq. (100) with respect to the angle O" to the third 
order and making a linear transformation like Eq. (97), one arrives at Eq. (59). We introduced in (59) 
several new quantities like 



Wj_ = 






l(3c3)4/3 



>0, 






/W_i_, 



7^ 



?(Wi)-" 



1 I T-L^* _ 1 / C" Oa „ () Da O Dq, 



FlM 



IM 



SttMKW ' 



(104) 
(105) 
(106) 



where Fim is the stability factor for long diameters, see Eq. (61), 
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Note, according to Eq. (106), the quantity c^ goes to near the separatrix {a -^ 1) like for the caustic 
case. This is the reason why the Maslov-Fedoryuk theory [39, 40, 41, 42] can be used for the transfor- 
mation of the integral over angle O" in Eqs. (100) into Eq. (59). 



C Jacobians for closed orbits with two reflection points 

The Jacobian J'^^^ defined by the derivative in Eq. (64) for closed orbits a like "co2" with two reflection 

points, J^^2 = {^y"/^^p) 2' '^^^ ^^ calculated by means of the caustic method [11]. The main idea of 
this method is to use a specific property of the trajectories in billiard system like elliptic cavity. They are 
the straight lines which tangent a curve called the elliptic or hyperbolic caustics between turning points. 
Our trajectory stability problem for the variations Sy" at a given S9p, see Fig. 3, is much simplified by 
reducing it to the calculation of the caustics semi- axes Uc, be and Uc -|- Suc, be + Sbc for closed orbit "co2" 
and its 66' deflection, respectively. For the case of closed non-periodic orbits "co2" the semi-axes ac and 
be and their variations are functions of the initial point (x, y) in contrast to the stability problem for the 
periodic orbits of Ref. [11]. The orbit-length invariant curve (confocal-to-boundary ellipse or hyperbola 
crossing the point {x,y), see Fig. 4) and its semi-axis variations play a similar role for the calculation of 
the "co2" stability factor jj}^2 with that of the boundary parameter for the periodic orbits in Ref. [11]. 
In this way this stability factor is obtained in the form 
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where x" is the x-coordinate of the final point O' (see Fig. 3), qo and qi are the tangents of the slope 
angle for the initial and final directions of particle motion along the orbit "co2" , 
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Here, the upper and lower signs stand for the hyperbolic and elliptic closed orbits, (xd, yd) and {xc2, yc2) 
are the first and last tangent-to-caustics points of the trajectory "co2" , 
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respectively, and 

Ac = blx' T aly\ Be = Talh\x, Cc - a\(hl - y"). 

The semi-axes Oc and he as functions of the initial point (x, y) for the hyperbolic or elliptic caustics for 
the orbit "co2" (sec Fig. 4) are given by 
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where ax and bx are the semi-axes for the confocal-to-boundary hyperbola or the ellipse crossing any 
current initial and final point {x, y) of the orbit "co2" inside the elliptic billiard. 
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and Z is the root of the cubic algebraic equation. 
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The factor T) in Eq. (108) is given by 
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We used here the invariancc of the Jacobian J^{x, y) against time reversal. 
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